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We use a multifractal formalism to study the effect of stochastic resonance in a noisy bistable 
system driven by various input signals. To characterize the response of a stochastic bistable system 
we introduce a new measure based on the calculation of a singularity spectrum for a return time 
sequence. We use wavelet transform modulus maxima method for the singularity spectrum com- 
putations. It is shown that the degree of multifractality defined as a width of singularity spectrum 
can be successfully used as a measure of complexity both in the case of periodic and aperiodic 
(stochastic or chaotic) input signals. We show that in the case of periodic driving force singularity 
spectrum can change its structure qualitatively becoming monofractal in the regime of stochastic 
synchronization. This fact allows us to consider the degree of multifractality as a new measure of 
stochastic synchronization also. Moreover, our calculations have shown that the effect of stochastic 
resonance can be catched by this measure even from a very short return time sequence. We use also 
the proposed approach to characterize the noise-enhanced dynamics of a coupled stochastic neurons 
model. 
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I. INTRODUCTION 

It is well known that noise is present inevitably in all 
real processes. The reckoning of its influence is very im- 
portant for a deeper understanding of the dynamics of 
real systems . Stochastic resonance (SR) discovered by 
Benzi, et al. and Nicolis et al. JlJ during the study of 
the Ice Ages is one of the bright examples of a nontrivial 
noise action on a nonlinear system. As a model of climate 
dynamics they proposed to consider a bistable system si- 
multaneously driven by noise and a periodic signal. It 
was shown that tuning the noise level, an enhancement 
of a bistable system's response to the periodic force be- 
comes possible. Beginning from the late 1980s, a wealth 
of theoretical and experimental papers followed, extend- 
ing the notion of SR (for extensive reviews, see ) and 
discovering new applications in different fields of sciences. 
There are a few different approaches to quantitative de- 
scription of stochastic resonance depending on the ampli- 
tude and character of an input signal which may be peri- 
odic, chaotic or even stochastic. Originally, the enhance- 
ment of a weak periodic input signal was characterized 
by the response amplitude at the frequency of periodic 
signal. Fauve and Heslot [^) and McNamara, Wiesen- 
feld, and Roy Q suggested to use the signal-to-noise ratio 
(SNR) as a quantitative measure of SR. Both quantities, 
the amplitude and the SNR, undergo a resonance-like 
curve as a function of noise intensity. Spectrum power 
amplification defined in ||J^] as the ratio of periodic com- 
ponents in output and input power spectrums demon- 
strates similar behavior taking a maximal value at an 
optimal noise level. 

Other measures based on the residence-time distribu- 
tion were introduced for description of SR in H^]. In 
this case the main object of considerations is the struc- 



ture of the mentioned distribution that contains a series 
of peaks at the odd multiples of the half-period of driv- 
ing. All of them go through maxima as a function of the 
noise strength || . Gammaitoni et al. introduced the area 
under the peak of the residence-time distribution at the 
half-period of driving as a measure of SR |l(]] . Some later, 
a fully systematic theory for the residence-time distribu- 
tion functions was developed by Choi, Fox and Jung [ pd| . 
They showed that to characterize correctly SR based on 
the residence-time distribution it is necessary to find the 
difference between the residence-time distribution in the 
presence of the modulation and the residence-time dis- 
tribution in the absence of the modulation at the half- 
period of the external force. In [|l2| the receiver operating 
characteristic was used for quantitative characterization 
of a response of coupled overdamped nonlinear dynamic 
elements driven by a weak sinusoidal signal embedded 
in Gaussian white noise. This approach was comple- 
mented and generalized in jl3| where SR was described in 
terms of maximization of information-theoretic distance 
measures between probability distributions of the output 
variable. 

Evidently, the measures mentioned above can be used 
in the case when the input signal has a clear distin- 
guishable peak in its power spectrum. To characterize 
a response of a noisy nonlinear system on an aperiodic 
driving force it is necessary to use other measures. In 
order to estimate the response of a noisy excitable (or 
bistable) system to a weak aperiodic signal Collins et al. 
|l4| introduced the input-output cross-correlation mea- 
sures and a measure of transinformation quantifying the 
rate of information transfer from stimulus to response. 
They showed, in particular, that the rate of information 
transfer between system output and input is optimized by 
noise and coined the term aperiodic stochastic resonance 



to describe this phenomenon. The coherence function 
was used in jl5| to characterize the response of a bistable 
system to a weak stochastic input. It was calculated an- 
alytically in the framework of the linear response theory 
(LRT) that has been successfully applied to SR and re- 
lated phenomena P{l(|. 

From the practical point of view, it is very impor- 
tant to have the measures calculated from a sequence 
of the time intervals characterizing the dynamics of an 
object under study. Two examples of such sequences, 
which may be most popular at present, are interbeat 
interval time series from cardiophysiology and neurons 
spike trains from neurodynamics. The statistical anal- 
ysis of heart-beat dynamics have shown that the use of 
approaches basing on wavelet or Hilbert transform (or 
their dual use) has a number of benefits in comparison 
with the traditional ones such as power spectrum and 
correlation analysis jL7|. It allows to analyze the in- 
formation stored in the Fourier phases of a signal un- 
der study which is crucial for determination of nonlin- 
ear characteristics. As was recently discovered by Ivanov 
et al. ]Tg| ] the human heartbeat dynamics possesses the 
multifractal properties. They used the wavelet-based ap- 
proach developed in []l9|-p5f to the analysis of complex 
non-stationary time series. It has been shown that a 
heartbeat sequence of a healthy subject has a multifractal 
scaling whereas the data from subjects with a pathologi- 
cal condition demonstrate a loss of mult ifr act ality. More- 
over, authors demonstrated an explicit relation between 
the nonlinear features (represented by the Fourier phase 
interactions) and the multifractality of healthy cardiac 
dynamics. The efficiency of their approach based on the 
time-frequency localization properties of wavelets allow- 
ing to analyze the non-stationarityies in time series. 

It is reasonable to try to use the same approach for 
the quantitative characterization of SR when a sequence 
under study is defeated by the concerted action of an 
external signal and of a random force. In this case, 
an input signal (periodic or aperiodic) is the source of 
non-stationarity in a response which can be analyzed by 
means of the wavelet-based algorithm mentioned above. 
From this point of view, it is naturally to expect that 
the use of the multifractal formalism will allow us to in- 
troduce a new universal measure quantifying SR for an 
arbitrary external signal. 

The main goal of the present study is the descrip- 
tion of SR from the multifractal analysis point of view. 
Basically, we treat as a model a bistable system simul- 
taneously driven by the white noise and an input sig- 
nal. To characterize the scaling properties of a return 
time sequence we use a spectrum of local Holder ex- 
ponents calculated by means of the wavelet transform 
modulus-maxima (WTMM) method. For this purpose, 
we first give some necessary definitions and illustrate the 
procedure of singularity spectrum calculation in Sec. II. 
Section III deals with the multifractal analysis of the 
stochastic bistable system response for different kinds of 
input signals. In this section, we introduce a new measure 



for quantitative description of SR and stochastic synchro- 
nization. We also test it ability to catch SR for different 
lengths of a return time sequence. In Sec. IV we apply 
mentioned approach to the study of aperiodic stochastic 
resonance and coherence resonance in an unidirection- 
ally coupled neurons model. In Sec. V we summarize 
our results and discuss the advantages of our measure in 
comparison with traditional measures. 



II. SINGULARITY SPECTRUM AND 
WAVELET-TRANSFORM MODULUS-MAXIMA 
METHOD 

It is well known that stochastic signals can be condi- 
tionally divided into two different classes. The first one 
includes the homogeneous signals characterizing by a sin- 
gle global Hurst exponent and having the same scaling 
properties at all time intervals. The second one includes 
the multifractal signals to describe the scaling properties 
of which it is necessary to use many local Hurst exponents 
(or Holder exponents) quantifying the local singular be- 
havior and local scaling in time series. According to the 
definition [ p5[ , the Holder exponent h(xo) of a function 
/ at the point xq is the greatest h so that / is Lipschitz 
at xq, i.e., there exists a constant C and a polynomial 
P n {x) of order n so that for all x in a neiborhood of x 
we have 

| f(x)-P n (x-x )\<C | x-x \ h . 

In fact, it measures the degree of irregularity of / at the 
point Xq. The singularity spectrum D{h) of the signal 
can be defined as the function that gives for a fixed h, 
the Hausdorff dimension of the set of points x where the 
exponent h(x) is equal to h. 

As was mentioned above, to determine the whole sin- 
gularity spectrum D{h) from an experimental signal it is 
necessary to use the approach based on the wavelet trans- 
form (WT), which permits an analysis both in physical 
space and in scale space. The WT of the function / is 
defined as 



Tih{b, a) = - 
a 



f(x) dx, 



(1) 



where ip is the analyzing wavelet, a € R + is a scale pa- 
rameter and b € R is a space parameter. The analyzing 
wavelet tp is generally chosen to be well localized in both 
space and frequency domain. A class of widely used real- 
valued analyzing wavelets which satisfies the above condi- 
tion is given by the Gaussian function and its derivatives. 
As was proven by Mallat and Hwang 21 the WT modu- 
lus maxima (local maxima of | T^{x, a) | at a given scale 
a) detect all the singularities of a signal under study. The 
skeleton from the modulus maxima lines contains all the 
information about the hierarchical distribution of singu- 
larities in the signal. The WTMM method consists in 
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taking advantage of the space-scale partitioning given by 
this skeleton to define a partition function which scales, 
in the limit a — ► + , in the following way [^,^3 24 : 



Z(q, a) 



E 



T^(xt(a),a) a T 



(2) 



where {xi(a)}i are the WT modulus maxima and q G R. 
According to the theorem proved in p4| . D{h), the singu- 
larity spectrum of the function /, is obtained by Legendre 
transformation of the function r(q) defined in (^|) 

D(h) = m\n(qh — r(q)). 
i 

The variables h and D(h) play the same role as the energy 
and entropy in the thermodynamics, whereas instead of 
the inverse of temperature and free energy we have q and 
r(q) |23| . |2"5|| . From a numerical point of view, it is more 
conveniently to calculate at first the scaling exponents: 

h(q) = lim — !- V* TU(g; ffj(a), a) In | T^(xi(a), a) |, 



{x z (a)} z 



and 



D{h(q)) = lim — !- T,/,(g; ^(o), a) lnT^g; £;(a), a), 

n — >n n n ^ — ^ 



0^0 In a 

{xi(a)h 

where f^,(q;x t (a),a) =| T^(xi,a) \ q /J2 Xi I ( x « > a ) I 9 - 
Further, we extract the set of Holder exponents and cor- 
responding singularity spectrum D{h) from log-log plots 
of h(q) and £>(%)) ||. 

As a simple example, we calculated r{q) and D(h) for 
the ordinary Brownian motion which is characterized by 
the single global Hurst's exponent H = 1/2. Figure 1 
demonstrates clearly the homogeneous scaling for the or- 
dinary Brownian motion. For more detailed information 
about calculation's procedure and additional references 
on freely distributed software see pq|. 
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FIG. 1. (a) - the path of a Brownian particle; (b) - the 
modulus-maxima skeleton of the random signal pictured in 
(a); (c) - the dependence r(q); (d) - the singularity spec- 
trum. The first derivative of the Gaussian function was used 
as the analyzing wavelet. 



III. MULTIFRACTAL APPROACH TO 
STOCHASTIC RESONANCE 



We treat for our study the overdamped bistable os- 
cillator which is governed in canonical units [Q] by the 
following stochastic differential equation (SDE): 



x = x — x 



2Q£{t)+y(t), 



(3) 



where is the white Gaussian noise with the correla- 
tion function (£(*)£(*')) = 5(t-t') and (£(*)) = 0, y(t) 
is an input signal. 

In the following subsections, we present the results of 
numerical simulations of Eq. (^) for different kinds of 
input signal. To characterize SR from the multifractal 
formalism point of view we will process the sequences 
of the return times to the one of the potential wells 0] 
normalized on a characteristic time scale of an external 
driving force. 



A. Periodic input signal 



Let us start from the more simple and well studied case 
of periodic external driving when y(t) — A sin(f2t + (f>). 
The amplitude of the input signal is assumed to be small, 
i.e., the signal alone cannot switch the system from one 
state to another in the absence of noise. For the low- 
frequency periodic modulation considered in this paper 
this means 



A < An = 



3V3' 



(4) 



As is well known, the study of SR can be conditionally 
separated on the two cases. The first one deals with the 
situation of a weak input signal when the amplitude of 
periodic driving is very small in comparison with a po- 
tential barrier and can be considered in the framework of 
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the LRT |T(j|n|. The second one is beyond the limits of 
LRT and corresponds to an amplitude of a subthreshold 
input signal comparable with a barrier. In the last case, 
the dynamics of a bistable system is characterized by a 
high degree of coherence between the switching process 
and input signal j27|]. It can be correctly described in 
terms of the phase synchronization theory both for pe- 
riodic | p8[ , chaotic |2{| and stochastic input signals fl30| . 
Moreover, as was shown in |3l| SR takes the form of 
the noise-induced order in this case. Such information- 
theoretical measures as the source entropy and the dy- 
namical entropy display a minimum at an optimal noise 
intensity. 
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FIG. 2. (a) D(h) in the case of the weak periodic input 
signal A — 0.081, Q = 0.004 for different values of noise 
intensity Q; (b) r(q) corresponding to the singularity spec- 
trums, D(h), in (a) for different values of noise intensity; 
the return time sequences normalized on the external force 
period and corresponding singularity spectrums in the case 
of A = 0.286, n = 0.004 for Q = 0.040 ((c) and (d)) and 
Q = 0.012 ((e) and (f)). The Gaussian function was used 
as the analyzing wavelet. 

In order to describe quantitatively these different situa- 
tions from the multifractal formalism point of view, we 
calculated singularity spectrums for different values of 
the driving amplitude taking as the signal under study 
a sequence of the return times normalized on the driv- 
ing period. Each sequence contained 10000 points. The 
results of calculations have shown the following. Singu- 



larity spectrum of the response has a bell shape form 
both in the presence and in the absence of the exter- 
nal periodic driving. That caused by a nonlinearity of 
the bistable system's response to the external driving 
force which manifests itself in the presence of the dif- 
ferent modes and in their interaction. As seen from 
Fig. 2, the tuning of noise level in the system of Eq. 
(||) leads to the changes in singularity spectrum both for 
weak and for strong enough driving signals. The width 
of the singularity spectrum takes its minimal value for 
an optimal noise intensity. It remains finite in the case 
of a weak periodic driving force, whereas in the regime of 
stochastic synchronization singularity spectrum qualita- 
tively changes its form shrinking to a single point. The 
return times fluctuate around the driving period in the 
regime of noise-enhanced phase locking and large bursts 
are seldom happen (Fig. 2 (c)), while for the values of 
noise intensities lying outside of synchronization region 
the respectively large fluctuations dominate (Fig. 2 (e)). 
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FIG. 3. The degree of multifractality, /3, vs. noise in- 
tensity for different values of the periodic force amplitude: 
A = (stars), A = 0.081 (triangles), A = 0.152 (dia- 
monds), A = 0.268 (circles). The frequency of external 
force has the same value as in Fig. [] The Gaussian function 
was used as the analyzing wavelet. 

It is natural to consider the degree of multifractal- 
ity (3 defined as the difference between the maximal and 
minimal Holder exponents belonging to the one and the 
same singularity spectrum as a new measure for SR and 
stochastic synchronization. As clearly seen from Fig. 3, 
the dependence of (3 on the noise intensity is character- 
ized by the presence of a minimum both for weak and for 
sufficiently strong input signals. In the regime of stochas- 
tic synchronization (3 equals to zero that corresponds to 
a single point singularity spectrum (see Fig. 2 (d)). In 
this case the scaling features of the return time sequence 
under study is characterized by the single scaling expo- 
nent that caused by the linearization of the response in 
the regime of switchings synchronization. Indeed, the 
interaction between different modes in response is sup- 
pressed in synchronous regime, because the switchings in 
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the system (|^) are in-phase with the input signal. The 
mode corresponding to the periodic input signal domi- 
nates and suppresses all others. 
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FIG. 5. The degree of multifractality vs. noise intensity 
calculated for indicated values of the length of the return 
time sequence. The amplitude and frequency of the driv- 
ing force are A = 0.12, Q = 0.004. 
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FIG. 4. Regions of stochastic synchronization con- 
structed by means of effective diffusion constant (dashed 
line) and of multifractality degree (solid line) . Driving fre- 
quency is fl = 0.004. The Gaussian function was used as 
the analyzing wavelet. 



In turn, it leads to simplification of the singularity spec- 
trum that reflects the absence of interactions between 
Fourier phases in response. It is necessary to empha- 
size here, that we understand stochastic synchronization 
as instantaneous matching of the input/output phases 
which is observed in a finite region on the parameter 
plane "noise intensity - amplitude of periodic force" . 
Traditionally, synchronization in noisy nonlinear oscil- 
latory systems is estimated quantitatively by means of 
effective diffusion constant which characterizes a veloc- 
ity of spreading of an initial phase difference distribution 
|p2| . Recently, this classical approach to synchroniza- 
tion was successfully used for quantitative description of 
the noise-enhanced phase coherence which takes place 
in stochastic bistable system driven by a subthreshold 
external signal p8|-|30||. The effective diffusion constant 
demonstrates a minimum decreasing up to a very small 
value in the region of synchronization. The above intro- 
duced wavelet-based measure demonstrates exactly the 
same behavior taking the zero value in the phase-locking 
regime that allows us to consider it as the measure of 
stochastic synchronization as well. 



Using (3 as the measure of stochastic synchronization 
it is possible to construct the region of synchronization 
on the parameter plane "noise intensity - amplitude of 
periodic force" (see Fig. 4). Inside of this region singu- 
larity spectrum of the response remains monofractal that 
once more time demonstrates the simplification of the re- 
sponse in the regime of synchronization. Synchronization 
region has a tongue-like form and nearly coincides with 
the similar one constructed by means of the effective dif- 
fusion constant used in Pq-pO||. It should be noted that 
the degree of multifractality does not relate directly with 
effective diffusion constant. The multifractal approach 
based on the analysis of the scaling features of the tem- 
poral sequence whereas the effective diffusion constant is 
the value characterizing a probability distribution of the 
input/output instantaneous phase difference. Our nu- 
merical studies also have shown that j3 demonstrates a 
monotonous dependence on the driving frequency. 



The length of the analyzed signals becomes the impor- 
tant parameter if they obtained in real experiments with 
live objects. Evidently, to estimate the possibility to use 
the above proposed measure in real situations we need 
to test its ability to catch SR for different lengths of the 
return time sequences. The results of our computations 
have been shown that the degree of multifractality takes 
the possibility to observe SR both for long and for suffi- 
ciently short return time sequences. As seen from Fig. 5, 
f3 calculated over the short return time sequences con- 
taining only 1000 points demonstrates nearly the same 
behavior as in the case of the long sequences. 
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FIG. 6. Spectral power amplification vs. noise inten- 
sity for different values of the driving amplitude. Driving 
frequency is f2 = 0.004. 

Thus, the width of singularity spectrum manifests itself 
as the general and universal characteristic for description 
of SR, because it works very well in a wide range of values 
of the driving amplitude and frequency. Its calculation 
allows us to get the information about the scaling prop- 
erties of the frequency fluctuations of a nonlinear system 
response and doesn't requires any information about in- 
put signal. Moreover, that new characteristic allows one 
to analyze effectively a response even in the case of short 
return time sequences that has an especial meaning for 
possible applications. 

The obtained results are in a good agreement with the 
behavior of the spectral power amplification tj firstly used 
for quantitative characterization of SR in pp] . ft demon- 
strates a maximum the absolute value of which is in- 
creased with the decrease of the input driving amplitude 
(see Fig. 6). The values of noise intensity maximizing 
r\ nearly coincide with the ones minimizing the degree 
of multifractality. The growth of the driving amplitude 
leads to the shift of minimal values of (3 to lower noise 
level as well as for spectral power amplification. 



B. Stochastic input signal 

From the practical application point of view, it is more 
interesting the situation when an input signal has a com- 
plex structure. SR for the input signals with fluctuating 
amplitude and phase was considered in |53|,[34| . In order 
to model the situation when the external signal is close 
to periodic one, but has a finite width of the spectral line 
Neiman and Schimansky-Geier [fi4| proposed to consider 
the harmonic noise as the input signal. Using the cumu- 
lant analysis and computer simulations they showed that 
the effect of SR takes place for harmonic noise as well 
and the width of the spectral line of the input signal at 
the output power spectrum can be decreased via SR. 

Harmonic noise y(t) i s defined by the following two- 
dimensional SDE |B3,P4|: 



where £(t) is the zero-mean Gaussian noise with 
(C(t)((t')} = 8{t — t'). It is necessary to note, that Gaus- 
sian noise Q(i) is statistically independent from the noise 
in (H). Equation (||) determines the two-dimensional 
Ornstein-Uhlenbek process y(t), s(t) with the power spec- 
trum 



Syy(uj) 



Ts 



lu 2 r 2 + (lu 2 - n 2 ) 2 ' 



(6) 



which for ft 2 > T 2 /4 has a peak at the frequency 
uj p = ^tt 2 - T 2 /2 with the width 



An 



uj 2 + Tlji - Jlj 2 - Tun, 



(7) 



where lj\ — y^ 2 " 
ments (y 2 ) = e/Vt 2 , 



- T 2 /4. The mean square displacc- 
(s 2 ) = e, (ys) = 0. The increase of 
the parameter F causes the widening of the spectral line 
p4j . We used harmonic noise as the input signal in (|^) 
to carry out the multifractal analysis of SR in the case 
of aperiodic driving force. As in the previous subsec- 
tion, we analyzed the return time sequences containing 
the same amount of points as before and normalized on 
the period T = 27r/0. To compare our results with those 
obtained previously, we choose the same values of param- 
eters for numerical simulations as in |34j ]. The results 
of our computations are presented in Fig. 7. The differ- 
ence between the maximal and minimal Holder exponents 
takes its minimal value for an optimal noise intensity as 
in the case of periodic driving. The obtained results are 
in good agreement with the results of 34 1. 
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FIG. 7. The degree of multifractality of the bistable sys- 
tem (^) response to the harmonic noise vs. noise intensity 
for the different values of the dissipation parameter in (g). 
Other parameters are Q = 0.1, e = 0.025. 
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The degree of multifractality behaves as the relative 
width of the output spectral line used in |j4] demonstrat- 
ing the minimum at a close value of the noise intensity. 
The decrease of F in (|J) leads to the regularization of 
the input signal that makes the minimum in f3(D) more 
pronounced. 



C. Chaotic input signal 

Now, let us take as the input signal in (||) the slowly 
varying subthreshold chaotic signal generated by the 
Lorenz system which is governed by the following ordi- 
nary differential equations: 



yi = 10 (2/2 - Vi)v, 

2/2 = (28 yi - j/2 - 2/12/3)^, 

2/3 = (2/12/2 ~ S/3y 3 )u, 



(8) 



where v is the small rationing constant slowing chaotic 
oscillations. Lorenz attractor existing in the phase space 
of this system has a thin multifractal structure [p5|. 
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FIG. 8. The input (a) and output (b) signals in the 
regime of stochastic synchronization for the case of chaotic 
driving signal generated by the Lorenz signal. The param- 
eters are k = 0.0188, Q = 0.03, v = 0.005. 
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FIG. 9. Singularity spectrums for different values of 
noise intensity Q and of the rationing constant: (a) 
k = 0.01, (b) k = 0.0188 
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FIG. 10. The degree of multifractality vs. noise inten- 
sity for the different values of chaotic signal amplitude. The 
width of the singularity spectrum of the input signal is rep- 
resented by the dashed line. The Gaussian function was 
used as the analyzing wavelet. 

The power spectrum of chaotic oscillations calculated for 
2/i,2-variables does not contain any sharp peaks in this 
case. The input signal has a form of random process 
of switchings between two metastable states (see Fig. 8 
(a)). It is reasonable to calculate the return time se- 
quences for the input and output signals and then to try 
to estimate the distortion of the signal in stochastic res- 
onator using singularity spectrum. We used the variable 
k ■ 2/1 (i) from the system of Eq. (||) as the input signal 
in our simulations, here k is a small positive rationing 
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constant. As follows from the results presented in Fig. 9 
(a), singularity spectrums of the input and output signals 
become more close to each other for some optimal noise 
level. Moreover, for some values of the parameter k the 
noise-enhanced phase coherence between chaotic signal 
and response is observed p9|] . In this case, the switch- 
ings in the input signal are in phase with the switchings 
in stochastic bistable system (see Fig. 8). Singularity 
spectrums of the input signal and of the response are co- 
incide (see Fig. 9 (b)) in the regime of phase locking that 
means the passing of chaotic signal through the stochas- 
tic resonator without any distortions. This is also illus- 
trated by the dependence of the singularity spectrum's 
width on the noise intensity for bistable system response 
that is presented in Fig. 10. For some optimal noise in- 
tensity, (3 takes its minimal value closing to the value 
corresponding to the width of the input signal's singu- 
larity spectrum. As was mentioned above, the input and 
output singularity spectrums coincide in the regime of 
stochastic synchronization that causes the coincidence of 
(3 with the multifractality degree of the external chaotic 
signal in some finite range of noise intensities. 

IV. MULTIFRACTAL ANALYSIS OF NEURON 
SPIKE TRAINS 

One of the reason for unremitting interest in SR is 
the possibility to model on its base different cooperative 
effects and the process of the information transfer in vari- 
ous biological sensory systems operating in natural noisy 
environment. At present, there are a lot of experimen- 
tal results showing that sensory neurons of different live 
organisms are able to demonstrate SR |pqfl]. In order 
to estimate the enhancement of a response, the signal- 
to-noise ratio or different cross-correlation measures are 
usually used |p^ , p6| . It is very interesting to use the 
above multifractal approach to analyze the spike trains 
generated by a stochastic neuronal model. 

We took as a model the Fitzhugh-Nagumo system |37j 
operating in excitable regime and driven by a mixture 
of the internal noise and a subthreshold stochastic spike 
train generated by another similar system detuned from 
the first one on a control parameter. The unidirectionally 
coupled neuron systems are described by the following 
stochastic differential equations: 

lix\ = xi - y - yi 

y\ = xi + ai + k x 2 + \/2Qi (9) 
fj,x 2 = x 2 - y - y 2 

y 2 = x 2 +a 2 + \/2Q 2 & (t) 

where a\ (Qi) and a 2 (Q 2 ) are, respectively, the control 
parameters (noise intensities) of the subsystems (x±,y±) 
and (x 2 ,y 2 ), £i and £2 are the statistically independent 



Gaussian white noise with the zero mean, A: is a small ra- 
tioning constant as before and ji -C 1 is a small parameter 
allowing one to separate all motions in the fast and slow 
ones. The values of control parameters and noise intensi- 
ties in subsystems are different and varied independently 
from each other. Thus, we can consider the system of 
Eq. (g) as a model of a single neuron embedded in a net- 
work and driven by both the internal noise and summed 
output of the neibouring neurons that can be modeled as 
a stochastic spike train. Interspike intervals (ISI) widely 
used in neuroscience as the typical neuron signals will 
play the role of signals under study in our consideration. 
The input stochastic spike train generated by the second 
neuron is characterized by a continuous singularity spec- 
trum having a finite width as well as the chaotic input 
signal from the Lorenz system. 

0.40 1 1 
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FIG. 11. Degree of multifractality of the first neuron re- 
sponse in Eq. (H) vs internal noise intensity Q\. Parame- 
ters: /i = 0.01, ai = 1.05, a 2 = 1.07, Q 2 = 0.02, k = 0.025. 
The width of the input signal singularity spectrum is la- 
beled by the dashed line. The Gaussian function was used 
as the analyzing wavelet. 

The calculated widths of the singularity spectrum for 
the ISI generated by the first neuron are shown in Fig. 11. 
It can be seen that the dependence of j3 on the inter- 
nal noise intensity is characterized by two different mini- 
mums corresponding to two different effects taking place 
in the system @ for small and large internal noise level, 
respectively. The first minimum appearing at the com- 
paratively small noise intensity corresponds to the effect 
of aperiodic stochastic resonance when a weak internal 
noise enhances the response of the neuron model opti- 
mizing the transmission of the input signal. As clearly 
seen from Fig. 11, the width of singularity spectrum cal- 
culated on the first neuron ISI is very close to the input 
one for some optimal values of noise level. Further in- 
crease of noise intensity makes singularity spectrum of 
response more narrow that can be considered as a mani- 
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festation of stochastic resonance without input signal pq ] 
called coherence resonance in Q. Indeed, in the case of 
sufficiently large noise intensity neuron already cannot 
distinguish the structure of the noisy input signal. It op- 
erates as an oscillator whose time scale is controlled by 
noise. For some optimal value of the internal noise in- 
tensity Q\ oscillations of the first neuron become close 
to the periodic one that is the essence of coherence res- 
onance [[39). At the moment when coherence resonance 
is observed f3 takes its second deeper minimum reflecting 
the noise-enhanced ordering of ISI. 



V. CONCLUSIONS 

We have studied the phenomenon of stochastic reso- 
nance in terms of the multifractal formalism revisited 
with wavelets p5[ . We observed that for some opti- 
mal noise intensity the degree of multifractality of the 
response, defined as a width of the singularity spectrum, 
takes its minimal value. Moreover, the qualitative change 
of its structure takes place in the regime of stochastic syn- 
chronization. In the region of the noise-enhanced phase 
locking it shrinks to the single point with zero Holder ex- 
ponent. We have shown that the width of the singularity 
spectrum calculated over the return time sequence can 
be effectively used as the measure characterizing the re- 
sponse of a noisy nonlinear system in a wide range of the 
driving amplitudes and frequencies. As follows from our 
numerical results this measure can be successfully used 
both for periodic and aperiodic (stochastic or chaotic) 
driving signals. Moreover, it has allowed us to estimate 
the degree of coherence for the unidirectionally coupled 
stochastic neurons model operating in excitable regime. 
By using the introduced measure, we successfully diag- 
nose both aperiodic stochastic resonance and coherence 
resonance which take place in the model under study for 
the small and large noise intensities, respectively. 

The proposed approach has a number of benefits in 
comparison with the traditionally used measures such as 
SNR, SPA, residence time distributions, coherence func- 
tion and others. These measures use the averaging pro- 
cedure for their calculation that leads to the loss of in- 
formation about nonlinear interaction between Fourier 
phases in response. This information is very important 
both for deeper understanding of the essence of SR and 
for more sensitive diagnostic of SR in full-scale experi- 
ments. The multifractal formalism based on wavelet cal- 
culations allows to study the temporal structure of the 
response. It catches all even weak non-stationarityies in a 
return times sequence under study that makes it a very 
powerful tool for diagnostic of SR and stochastic syn- 
chronization. The introduced measure demonstrates the 
behavior which in a very good agreement with the be- 
havior of traditional quantitative characteristics of SR. 
It is universal in relation to the kind of the input signal 
and able to catch noise-induced effects even from very 



short time series. Later has the especial importance for 
the analysis of real signals. 

The presented approach to the study of scaling features 
of motion in stochastic systems may be very fruitful also 
in the case of the Brownian motion in periodic potential 
under the action of random forces. It will be the task of 
our future investigations. 
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FIGURE CAPTIONS: 



Fig. 1. (a) - the path of a Brownian particle; (b) - the modulus-maxima skeleton of the random signal pictured in 

(a) ; (c) - the dependence r(q); (d) - the singularity spectrum. The first derivative of the Gaussian function was used 
as the analyzing wavelet. 

Fig. 2. (a) D(h) in the case of the weak periodic input signal A = 0.081,57 = 0.004 for different values of noise 
intensity Q; (b) r(q) corresponding to the singularity spectrums, D(h), in (a) for different values of noise intensity; 
the sequences of the return times normalized on the external force period and corresponding singularity spectrums in 
the case of A = 0.286,0 = 0.004 for Q = 0.040 ((c) and (d)) and Q = 0.012 ((e) and (f)); The Gaussian function 
was used as the analyzing wavelet. 

Fig. 3. Degree of multifractality vs. noise intensity for different values of the periodic force amplitude: A = 
(stars), A = 0.081 (triangles), A = 0.152 (diamonds), A = 0.268 (circles). The frequency of external force has the 
same value as in Fig. 2. The Gaussian function was used as the analyzing wavelet. 

Fig. 4. Regions of stochastic synchronization constructed by means of effective diffusion constant (dashed line) 
and of multifractality degree (solid line). Driving frequency is O = 0.004. The Gaussian function was used as the 
analyzing wavelet. 

Fig. 5. Degree of multifractality vs. noise intensity calculated for indicated values of the length of the return time 
sequence. The amplitude and frequency of the driving force are A = 0.12, 57 = 0.004. 

Fig. 6. Spectral power amplification vs. noise intensity for different values of the driving amplitude. Driving 
frequency is 57 — 0.004. 

Fig. 7. Degree of multifractality of the bistable system (||) response to the harmonic noise vs. noise intensity for 
the different values of the dissipation parameter in (g|). Other parameters are 57 = 0.1, e = 0.025. 

Fig. 8. The input (a) and output (b) signals in the regime of stochastic synchronization for the case of chaotic 
driving signal generated by the Lorenz signal. The parameters are k = 0.0188, Q = 0.03, v — 0.005. 

Fig. 9. Singularity spectrums for different values of noise intensity Q and of the rationing constant: (a) k = 0.01, 

(b) k — 0.0188. The Gaussian function was used as the analyzing wavelet. 

Fig. 10. Degree of multifractality vs. noise intensity for the different values of chaotic signal amplitude. The width 
of the singularity spectrum of the input signal is represented by the dashed line. The Gaussian function was used as 
the analyzing wavelet. 

Fig. 11. Degree of multifractality of the first neuron response in Eq. (^) vs internal noise intensity Q\. Parameters: 
/i = 0.01, a\ = 1.05, a2 = 1.07, Q2 — 0.02, k — 0.025. The width of the input signal singularity spectrum is labeled 
by the dashed line. The Gaussian function was used as the analyzing wavelet. 
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